import math
from sympy import *

h = float(input('开挖深度(m): '))
gamma = float(input('土层重度(kN/m^3): '))
phi = float(input('内摩擦角(°): '))
c = float(input('粘聚力(kPa): '))
q0 = float(input('地面超载(kPa): '))


rada = math.radians(45-phi/2)
radp = math.radians(45+phi/2)
k_a = math.tan(rada)**2
k_p = math.tan(radp)**2

pab = k_a*(q0+gamma*h) - 2*c*k_a**0.5
ppb = 0 + 2*c*k_p**0.5
pb1 = pab - ppb

u = (k_a*(q0+gamma*h) - 2*c*(k_p**0.5+k_a**0.5))/(gamma*(k_p-k_a))

z0 = 2*c/(gamma*k_a**0.5)-q0/gamma
ea1 = 0.5*pab*(h-z0)
ha1 = z0+2/3*(h-z0)
ea2 = 0.5*pb1*u
ha2 = h + 1/3*u
ea = ea1+ea2
ha = (ea1*ha1+ea2*ha2)/ea
x = Symbol('x')
t = solve(x**3 - 6*ea*x/gamma/(k_p-k_a) - 6*(h+u-ha)*ea/gamma/(k_p-k_a), x)
t = float(t[0])

xm = (2*ea/(k_p-k_a)/gamma)**0.5
m_max = (h+u+xm-ha)*ea-gamma*(k_p-k_a)*xm**3/h
print(xm,m_max)

print('悬臂式地下连续墙的有效嵌固深度为 %f m，最大弯矩为 %f kN·m' %(t, m_max))
